Input Data

DateTime Junction Vehicles ID
2015-11-01 00:00:00 1 15 20151101001
2015-11-01 01:00:00 1 13 20151101011
2015-11-01 02:00:00 1 10 20151101021
2015-11-01 03:00:00 1 7 20151101031
2015-11-01 04:00:00 1 9 20151101041
2015-11-01 05:00:00 1 6 20151101051
2015-11-01 06:00:00 1 9 20151101061
2015-11-01 07:00:00 1 8 20151101071
2015-11-01 08:00:00 1 11 20151101081
2015-11-01 09:00:00 1 12 20151101091
DateTime Junction Vehicles ID
48111 2017-06-30 14:00:00 4 10 20170630144
48112 2017-06-30 15:00:00 4 14 20170630154
48113 2017-06-30 16:00:00 4 16 20170630164
48114 2017-06-30 17:00:00 4 16 20170630174
48115 2017-06-30 18:00:00 4 17 20170630184
48116 2017-06-30 19:00:00 4 11 20170630194
48117 2017-06-30 20:00:00 4 30 20170630204
48118 2017-06-30 21:00:00 4 16 20170630214
48119 2017-06-30 22:00:00 4 22 20170630224
48120 2017-06-30 23:00:00 4 12 20170630234
## next
knitr::kable(rbind(head(cars),tail(cars)))
cat("next")
Date Hour junction1 junction2 junction3
1 2015-11-01 0 15 6 9
2 2015-11-01 1 13 6 7
3 2015-11-01 2 10 5 5
4 2015-11-01 3 7 6 1
5 2015-11-01 4 9 7 2
6 2015-11-01 5 6 2 2
14587 2017-06-30 18 95 34 38
14588 2017-06-30 19 105 34 33
14589 2017-06-30 20 96 35 31
14590 2017-06-30 21 90 31 28
14591 2017-06-30 22 84 29 26
14592 2017-06-30 23 78 27 39
## next
# knitr::kable(head(traffic.data))
# knitr::kable(head(traffic.data[traffic.data$Junction == 2,]))
# knitr::kable(head(traffic.data[traffic.data$Junction == 3,]))
 knitr::kable(head(traffic.data[traffic.data$Junction == 4,1:3]))

knitr::kable(head(cbind(traffic.data[traffic.data$Junction == 1,1:3][1:10,],
       traffic.data[traffic.data$Junction == 2,2:3][1:10,],
       traffic.data[traffic.data$Junction == 3,2:3][1:10,])))

hold = cbind(traffic.data[traffic.data$Junction == 1,1:3],
             traffic.data[traffic.data$Junction == 2,2:3],
             traffic.data[traffic.data$Junction == 3,2:3])

cat("next")
DateTime Junction Vehicles
43777 2017-01-01 00:00:00 4 3
43778 2017-01-01 01:00:00 4 1
43779 2017-01-01 02:00:00 4 4
43780 2017-01-01 03:00:00 4 4
43781 2017-01-01 04:00:00 4 2
43782 2017-01-01 05:00:00 4 1
DateTime Junction Vehicles Junction Vehicles Junction Vehicles
2015-11-01 00:00:00 1 15 2 6 3 9
2015-11-01 01:00:00 1 13 2 6 3 7
2015-11-01 02:00:00 1 10 2 5 3 5
2015-11-01 03:00:00 1 7 2 6 3 1
2015-11-01 04:00:00 1 9 2 7 3 2
2015-11-01 05:00:00 1 6 2 2 3 2
## next

create daily cars

dates = unique(cars$Date)
cars.day = data.frame(Date = dates, Junction1 = rep(NA,length(dates)), Junction2 = rep(NA,length(dates)), Junction3 = rep(NA,length(dates)))

for(i in 1:length(dates)){
  cars.day$Junction1[i] = sum(cars$junction1[cars$Date == dates[i]])
  cars.day$Junction2[i] = sum(cars$junction2[cars$Date == dates[i]])
  cars.day$Junction3[i] = sum(cars$junction3[cars$Date == dates[i]])
}

cars.day$Junction1 = ts(cars.day$Junction1, frequency = 7)
cars.day$Junction2 = ts(cars.day$Junction2, frequency = 7)
cars.day$Junction3 = ts(cars.day$Junction3, frequency = 7)

knitr::kable(rbind(head(cars.day),tail(cars.day)))

plot(cars.day$Junction1, type = "n", ylim = c(0,2250), main = "daily cars at junctions", xlab = "Time (day)", ylab = "cars")
grid()
lines(cars.day$Junction1, col = rainbow(3)[1])
lines(cars.day$Junction2, col = rainbow(3)[2])
lines(cars.day$Junction3, col = rainbow(3)[3])
legend("topleft",legend = c("Junction 1", "Junction 2", "Junction 3"), col = rainbow(3),lty = c(1,1,1))

plot(cars.day$Junction1, type = "n", ylim = c(0,2250), main = "daily cars at junctions 1 & 2", xlab = "Time (day)", ylab = "cars")
grid()
lines(cars.day$Junction1, col = rainbow(3)[1])
lines(cars.day$Junction2, col = rainbow(3)[2])
legend("topleft",legend = c("Junction 1", "Junction 2"), col = rainbow(3)[1:2],lty = c(1,1))
Date Junction1 Junction2 Junction3
1 2015-11-01 327 133 136
2 2015-11-02 546 197 166
3 2015-11-03 544 217 150
4 2015-11-04 498 199 121
5 2015-11-05 464 200 106
6 2015-11-06 446 199 101
603 2017-06-25 1184 367 335
604 2017-06-26 1774 587 414
605 2017-06-27 2187 730 567
606 2017-06-28 2080 720 503
607 2017-06-29 2086 696 497
608 2017-06-30 1883 655 553

acf(cars.day$Junction1, lag.max = 30, main = "Junction 1 ACF")
pacf(cars.day$Junction1, lag.max = 30, main = "Junction 1 PACF")

acf(cars.day$Junction2, lag.max = 30, main = "Junction 2 ACF")
pacf(cars.day$Junction2, lag.max = 30, main = "Junction 2 PACF")

acf(cars.day$Junction3, lag.max = 30, main = "Junction 3 ACF")
pacf(cars.day$Junction3, lag.max = 30, main = "Junction 3 PACF")

create traning data

# remove last month data

nday = nrow(cars.day)

cars.day.train = cars.day[-((nday-29):nday),] # remove last month

cars.day.train$Junction1 = ts(cars.day.train$Junction1,frequency = 7)
cars.day.train$Junction2 = ts(cars.day.train$Junction2,frequency = 7)
cars.day.train$Junction3 = ts(cars.day.train$Junction3,frequency = 7)

auto arima

Junction 1

(cars.day.fit1 = auto.arima(cars.day.train$Junction1,
                           max.p = 5,
                        max.q = 5,
                        max.P = 5,
                        max.Q = 5,
                        max.order = 10,
                        max.d = 3,
                        max.D = 3,
                        start.p = 1,
                        start.q = 1,
                        start.P = 1,
                        start.Q = 1,
                        nmodels = 100,
                        approximation = T))


#quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d = 1, q = 2))
## Series: cars.day.train$Junction1 
## ARIMA(1,0,1)(0,1,1)[7] with drift 
## 
## Coefficients:
##          ar1      ma1     sma1   drift
##       0.8128  -0.3473  -0.8211  2.1295
## s.e.  0.0386   0.0629   0.0243  0.3114
## 
## sigma^2 = 6325:  log likelihood = -3310.96
## AIC=6631.91   AICc=6632.02   BIC=6653.65

Junction 2

(cars.day.fit2 = auto.arima(cars.day.train$Junction2,
                           max.p = 3,
                        max.q = 3,
                        max.P = 3,
                        max.Q = 3,
                        max.order = 5,
                        max.d = 1,
                        max.D = 1,
                        start.p = 1,
                        start.q = 1,
                        start.P = 1,
                        start.Q = 1,
                        nmodels = 100,
                        approximation = T))

quiet(astsa::sarima(cars.day.train$Junction2, p = 3, d = 1, q = 3))
## Series: cars.day.train$Junction2 
## ARIMA(2,0,1)(0,1,1)[7] with drift 
## 
## Coefficients:
##          ar1      ar2      ma1     sma1   drift
##       1.4178  -0.4346  -0.8493  -0.7995  0.7271
## s.e.  0.0652   0.0593   0.0444   0.0256  0.2996
## 
## sigma^2 = 778.1:  log likelihood = -2711.57
## AIC=5435.15   AICc=5435.3   BIC=5461.23

Junction 3

(cars.day.fit3 = auto.arima(cars.day.train$Junction3,
                           max.p = 3,
                        max.q = 3,
                        max.P = 3,
                        max.Q = 3,
                        max.order = 5,
                        max.d = 1,
                        max.D = 1,
                        start.p = 2,
                        start.q = 2,
                        start.P = 2,
                        start.Q = 2,
                        nmodels = 50,
                        approximation = T))

quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d = 1, q = 1))
## Series: cars.day.train$Junction3 
## ARIMA(1,1,1)(0,0,1)[7] 
## 
## Coefficients:
##          ar1      ma1    sma1
##       0.6969  -0.9755  0.0675
## s.e.  0.0354   0.0125  0.0420
## 
## sigma^2 = 8704:  log likelihood = -3435.03
## AIC=6878.07   AICc=6878.14   BIC=6895.5

find seasonal model

Junction 1

cat("model 1")
(model1.1 = quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d = 1, q = 2, P = 1, D = 0, Q = 0, S = 7)))
#AIC 6953.23

# model1.forecast2 = forecast(model1.1, h = 30)
# rmse1.1 = rmse(model1.forecast1$mean, data.1.full$J1[577:606])
# 114.3474


cat("model 2")
(model1.2 = quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d = 1, q = 2, P = 0, D = 0, Q = 1, S = 7)))
#AIC 7297.4

cat("model 3")
(model1.3 = quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d= 0, q = 1, P = 1, D = 1, Q = 0, S = 7)))
#AIC 6744.5

cat("model 4")
(model1.4 = quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d= 1, q = 0, P = 2, D = 1, Q = 0, S = 7)))
#AIC 6749.55

cat("model 5")
(model1.5 = quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d= 1, q = 0, P = 2, D = 0, Q = 0, S = 7)))
#AIC 6897.51


data.1 = data.frame(J1 = cars.day.train$Junction1[-c(1,2)],
                    J2L1 = cars.day.train$Junction2[-c(1,578)],
                    J2L2 = cars.day.train$Junction2[-c(577,578)],
                    J3L1 = cars.day.train$Junction3[-c(1,578)],
                    J3L2 = cars.day.train$Junction3[-c(577,578)])

data.1.full = data.frame(J1 = cars.day$Junction1[-c(1,2)],
                         J2L1 = cars.day$Junction2[-c(1,608)],
                         J2L2 = cars.day$Junction2[-c(607,608)],
                         J3L1 = cars.day$Junction3[-c(1,608)],
                         J3L2 = cars.day$Junction3[-c(607,608)])

cat("model 6")


# sarima(data.1$J1, p = 1, d= 1, q = 0, P = 2, D = 0, Q = 0, S = 7, 
#        xreg = cbind(data.1$J2L1,data.1$J3L1))

(model1.6 = quiet(astsa::sarima(data.1$J1, 
                              p = 1, 
                              d= 1, 
                              q = 0, 
                              P = 2, 
                              D = 0, 
                              Q = 0, 
                              S = 7, 
                              xreg = data.1[,-1])))

#AIC 6871.87

# model1.forecast1 = forecast(model1.6, h = 30)
# rmse1.1 = rmse(model1.forecast1$mean, data.1.full$J1[577:606])
# 114.3474

cat("model 7")
(model1.7 = quiet(astsa::sarima(data.1$J1, 
                              p = 1, 
                              d= 1, 
                              q = 1, 
                              P = 2, 
                              D = 0, 
                              Q = 0, 
                              S = 7, 
                              xreg = data.1[,-1])))

#AIC 6821.64

cat("model 8")
(model1.8 = quiet(astsa::sarima(data.1$J1, 
                              p = 2, 
                              d= 1, 
                              q = 1, 
                              P = 2, 
                              D = 0, 
                              Q = 0, 
                              S = 7, 
                              xreg = data.1[,-1])))

#AIC 6802.72

cat("model 9")
(model1.9 = quiet(astsa::sarima(data.1$J1, 
                              p = 2, 
                              d= 0, 
                              q = 1, 
                              P = 2, 
                              D = 1, 
                              Q = 0, 
                              S = 7, 
                              xreg = data.1[,-1])))

#AIC 6675.71

cat("model 10")
(model1.10 = quiet(astsa::sarima(data.1$J1, 
                              p = 2, 
                              d= 1, 
                              q = 1, 
                              P = 2, 
                              D = 1, 
                              Q = 0, 
                              S = 7, 
                              xreg = data.1[,-1])))

#AIC 6672.2

cat("model 11")
(model1.11 = quiet(astsa::sarima(data.1$J1, 
                              p = 3, 
                              d= 1, 
                              q = 1, 
                              P = 2, 
                              D = 1, 
                              Q = 0, 
                              S = 7, 
                              xreg = data.1[,-1])))

#AIC 6666.58

cat("model 12")
(model1.12 = quiet(astsa::sarima(data.1$J1, 
                              p = 2, 
                              d= 1, 
                              q = 3, 
                              P = 2, 
                              D = 1, 
                              Q = 0, 
                              S = 7, 
                              xreg = data.1[,-1])))

#AIC 6666.59

cat("model 13\n\n")
(model1.13<-Arima(data.1$J1, 
          order=c(1, 0, 1), 
          seasonal=list(order=c(0, 1, 1), period=7), 
          include.drift=T,
          xreg = as.matrix(data.1[,-1])))

model1.forecast1 = forecast(model1.13, h = 30, xreg = as.matrix(data.1.full[577:606,-1]))
rmse1.1 = rmse(model1.forecast1$mean, data.1.full$J1[577:606])
# 109.0436



#AIC 6610.41

cat("model 14n\n")
(model1.14<-Arima(data.1$J1, 
          order=c(1, 0, 1), 
          seasonal=list(order=c(0, 1, 1), period=7), 
          include.drift=T,
          xreg = as.matrix(data.1[,-c(1,3,4,5)])))

#AIC 6604.65


#plot(data.1$J2L1, data.1$J1)
## model 1$fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     xreg = constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1      ma1     ma2    sar1  constant
##       0.6072  -1.1318  0.1318  0.8665    2.1810
## s.e.  0.0632   0.0774  0.0773  0.0219    0.3239
## 
## sigma^2 estimated as 9633:  log likelihood = -3470.62,  aic = 6953.23
## 
## $degrees_of_freedom
## [1] 572
## 
## $ttable
##          Estimate     SE  t.value p.value
## ar1        0.6072 0.0632   9.6104  0.0000
## ma1       -1.1318 0.0774 -14.6256  0.0000
## ma2        0.1318 0.0773   1.7054  0.0887
## sar1       0.8665 0.0219  39.6120  0.0000
## constant   2.1810 0.3239   6.7330  0.0000
## 
## $AIC
## [1] 12.05066
## 
## $AICc
## [1] 12.05084
## 
## $BIC
## [1] 12.09597
## 
## model 2$fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     xreg = constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1      ma1      ma2    sma1  constant
##       0.1312  -0.5222  -0.4778  0.5960    2.1702
## s.e.  0.0651   0.0551   0.0549  0.0294    0.0884
## 
## sigma^2 estimated as 17586:  log likelihood = -3642.7,  aic = 7297.4
## 
## $degrees_of_freedom
## [1] 572
## 
## $ttable
##          Estimate     SE t.value p.value
## ar1        0.1312 0.0651  2.0165  0.0442
## ma1       -0.5222 0.0551 -9.4789  0.0000
## ma2       -0.4778 0.0549 -8.6971  0.0000
## sma1       0.5960 0.0294 20.2882  0.0000
## constant   2.1702 0.0884 24.5467  0.0000
## 
## $AIC
## [1] 12.64714
## 
## $AICc
## [1] 12.64732
## 
## $BIC
## [1] 12.69245
## 
## model 3$fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     xreg = constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1      ma1     sar1  constant
##       0.7900  -0.3288  -0.5089    2.2030
## s.e.  0.0393   0.0579   0.0367    1.1115
## 
## sigma^2 estimated as 7724:  log likelihood = -3367.25,  aic = 6744.5
## 
## $degrees_of_freedom
## [1] 567
## 
## $ttable
##          Estimate     SE  t.value p.value
## ar1        0.7900 0.0393  20.0824   0.000
## ma1       -0.3288 0.0579  -5.6820   0.000
## sar1      -0.5089 0.0367 -13.8661   0.000
## constant   2.2030 1.1115   1.9820   0.048
## 
## $AIC
## [1] 11.81174
## 
## $AICc
## [1] 11.81186
## 
## $BIC
## [1] 11.8498
## 
## model 4$fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##           ar1     sar1     sar2
##       -0.3676  -0.6937  -0.3084
## s.e.   0.0391   0.0400   0.0401
## 
## sigma^2 estimated as 7964:  log likelihood = -3370.77,  aic = 6749.55
## 
## $degrees_of_freedom
## [1] 567
## 
## $ttable
##      Estimate     SE  t.value p.value
## ar1   -0.3676 0.0391  -9.3900       0
## sar1  -0.6937 0.0400 -17.3483       0
## sar2  -0.3084 0.0401  -7.6968       0
## 
## $AIC
## [1] 11.84131
## 
## $AICc
## [1] 11.84139
## 
## $BIC
## [1] 11.87181
## 
## model 5$fit
## 
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S), 
##     xreg = constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##         REPORT = 1, reltol = tol))
## 
## Coefficients:
##           ar1    sar1    sar2  constant
##       -0.3652  0.4476  0.5056    2.4760
## s.e.   0.0391  0.0357  0.0365   45.8665
## 
## sigma^2 estimated as 8697:  log likelihood = -3443.76,  aic = 6897.51
## 
## $degrees_of_freedom
## [1] 573
## 
## $ttable
##          Estimate      SE t.value p.value
## ar1       -0.3652  0.0391 -9.3386   0.000
## sar1       0.4476  0.0357 12.5358   0.000
## sar2       0.5056  0.0365 13.8638   0.000
## constant   2.4760 45.8665  0.0540   0.957
## 
## $AIC
## [1] 11.95409
## 
## $AICc
## [1] 11.95421
## 
## $BIC
## [1] 11.99186
## 
## model 6$fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##     REPORT = 1, reltol = tol))
## 
## Coefficients:
##           ar1    sar1    sar2    J2L1     J2L2     J3L1     J3L2
##       -0.3823  0.4515  0.4979  0.2207  -0.2148  -0.0083  -0.0106
## s.e.   0.0404  0.0359  0.0368  0.1195   0.1171   0.0331   0.0330
## 
## sigma^2 estimated as 8587:  log likelihood = -3427.93,  aic = 6871.87
## 
## $degrees_of_freedom
## [1] 568
## 
## $ttable
##      Estimate     SE t.value p.value
## ar1   -0.3823 0.0404 -9.4714  0.0000
## sar1   0.4515 0.0359 12.5916  0.0000
## sar2   0.4979 0.0368 13.5462  0.0000
## J2L1   0.2207 0.1195  1.8458  0.0654
## J2L2  -0.2148 0.1171 -1.8336  0.0672
## J3L1  -0.0083 0.0331 -0.2500  0.8026
## J3L2  -0.0106 0.0330 -0.3221  0.7475
## 
## $AIC
## [1] 11.95108
## 
## $AICc
## [1] 11.95142
## 
## $BIC
## [1] 12.01166
## 
## model 7$fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##     REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1      ma1    sar1    sar2    J2L1    J2L2     J3L1   J3L2
##       0.5188  -0.9898  0.4956  0.4412  0.3452  0.1263  -0.0006  0.027
## s.e.  0.0429   0.0048  0.0375  0.0385  0.1358  0.1304   0.0341  0.034
## 
## sigma^2 estimated as 7880:  log likelihood = -3401.82,  aic = 6821.64
## 
## $degrees_of_freedom
## [1] 567
## 
## $ttable
##      Estimate     SE   t.value p.value
## ar1    0.5188 0.0429   12.0839  0.0000
## ma1   -0.9898 0.0048 -205.2275  0.0000
## sar1   0.4956 0.0375   13.2063  0.0000
## sar2   0.4412 0.0385   11.4501  0.0000
## J2L1   0.3452 0.1358    2.5421  0.0113
## J2L2   0.1263 0.1304    0.9688  0.3331
## J3L1  -0.0006 0.0341   -0.0168  0.9866
## J3L2   0.0270 0.0340    0.7948  0.4270
## 
## $AIC
## [1] 11.86373
## 
## $AICc
## [1] 11.86417
## 
## $BIC
## [1] 11.93188
## 
## model 8$fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##     REPORT = 1, reltol = tol))
## 
## Coefficients:
##         ar1     ar2      ma1    sar1    sar2    J2L1     J2L2    J3L1    J3L2
##       0.433  0.2043  -0.9928  0.4774  0.4735  0.3400  -0.0389  0.0046  0.0085
## s.e.  0.044  0.0438   0.0038  0.0369  0.0380  0.1271   0.1269  0.0336  0.0336
## 
## sigma^2 estimated as 7582:  log likelihood = -3391.36,  aic = 6802.72
## 
## $degrees_of_freedom
## [1] 566
## 
## $ttable
##      Estimate     SE   t.value p.value
## ar1    0.4330 0.0440    9.8509  0.0000
## ar2    0.2043 0.0438    4.6667  0.0000
## ma1   -0.9928 0.0038 -259.0924  0.0000
## sar1   0.4774 0.0369   12.9338  0.0000
## sar2   0.4735 0.0380   12.4478  0.0000
## J2L1   0.3400 0.1271    2.6753  0.0077
## J2L2  -0.0389 0.1269   -0.3066  0.7593
## J3L1   0.0046 0.0336    0.1376  0.8906
## J3L2   0.0085 0.0336    0.2519  0.8012
## 
## $AIC
## [1] 11.83083
## 
## $AICc
## [1] 11.83138
## 
## $BIC
## [1] 11.90655
## 
## model 9$fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##     REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1     ar2      ma1     sar1     sar2    J2L1    J2L2     J3L1
##       0.8444  0.0069  -0.4156  -0.6663  -0.2948  0.3507  0.0019  -0.0013
## s.e.  0.0947  0.0729   0.0947   0.0416   0.0407  0.1275  0.1260   0.0336
##         J3L2
##       0.0004
## s.e.  0.0335
## 
## sigma^2 estimated as 6996:  log likelihood = -3327.86,  aic = 6675.71
## 
## $degrees_of_freedom
## [1] 560
## 
## $ttable
##      Estimate     SE  t.value p.value
## ar1    0.8444 0.0947   8.9185  0.0000
## ar2    0.0069 0.0729   0.0941  0.9251
## ma1   -0.4156 0.0947  -4.3886  0.0000
## sar1  -0.6663 0.0416 -16.0083  0.0000
## sar2  -0.2948 0.0407  -7.2411  0.0000
## J2L1   0.3507 0.1275   2.7494  0.0062
## J2L2   0.0019 0.1260   0.0153  0.9878
## J3L1  -0.0013 0.0336  -0.0394  0.9686
## J3L2   0.0004 0.0335   0.0110  0.9912
## 
## $AIC
## [1] 11.73236
## 
## $AICc
## [1] 11.73293
## 
## $BIC
## [1] 11.8087
## 
## model 10$fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##     REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1     ar2      ma1     sar1     sar2    J2L1     J2L2     J3L1
##       0.4589  0.2249  -1.0000  -0.6481  -0.2838  0.2827  -0.0611  -0.0059
## s.e.  0.0437  0.0431   0.0047   0.0414   0.0409  0.1256   0.1253   0.0337
##          J3L2
##       -0.0024
## s.e.   0.0337
## 
## sigma^2 estimated as 7030:  log likelihood = -3326.1,  aic = 6672.2
## 
## $degrees_of_freedom
## [1] 559
## 
## $ttable
##      Estimate     SE   t.value p.value
## ar1    0.4589 0.0437   10.5116  0.0000
## ar2    0.2249 0.0431    5.2188  0.0000
## ma1   -1.0000 0.0047 -211.4468  0.0000
## sar1  -0.6481 0.0414  -15.6632  0.0000
## sar2  -0.2838 0.0409   -6.9340  0.0000
## J2L1   0.2827 0.1256    2.2502  0.0248
## J2L2  -0.0611 0.1253   -0.4881  0.6256
## J3L1  -0.0059 0.0337   -0.1757  0.8606
## J3L2  -0.0024 0.0337   -0.0718  0.9428
## 
## $AIC
## [1] 11.74683
## 
## $AICc
## [1] 11.74739
## 
## $BIC
## [1] 11.82327
## 
## model 11$fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##     REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1     ar2     ar3      ma1     sar1     sar2    J2L1     J2L2
##       0.4262  0.1714  0.1190  -1.0000  -0.6629  -0.2953  0.3334  -0.0248
## s.e.  0.0446  0.0471  0.0429   0.0047   0.0414   0.0409  0.1255   0.1246
##         J3L1     J3L2
##       0.0007  -0.0023
## s.e.  0.0336   0.0334
## 
## sigma^2 estimated as 6938:  log likelihood = -3322.29,  aic = 6666.58
## 
## $degrees_of_freedom
## [1] 558
## 
## $ttable
##      Estimate     SE   t.value p.value
## ar1    0.4262 0.0446    9.5476  0.0000
## ar2    0.1714 0.0471    3.6410  0.0003
## ar3    0.1190 0.0429    2.7756  0.0057
## ma1   -1.0000 0.0047 -214.3363  0.0000
## sar1  -0.6629 0.0414  -16.0078  0.0000
## sar2  -0.2953 0.0409   -7.2106  0.0000
## J2L1   0.3334 0.1255    2.6576  0.0081
## J2L2  -0.0248 0.1246   -0.1994  0.8420
## J3L1   0.0007 0.0336    0.0218  0.9826
## J3L2  -0.0023 0.0334   -0.0686  0.9454
## 
## $AIC
## [1] 11.73693
## 
## $AICc
## [1] 11.73762
## 
## $BIC
## [1] 11.82102
## 
## model 12$fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc, 
##     REPORT = 1, reltol = tol))
## 
## Coefficients:
##           ar1     ar2      ma1      ma2     ma3     sar1     sar2    J2L1
##       -0.1167  0.8114  -0.4463  -0.9757  0.4221  -0.6771  -0.2916  0.3528
## s.e.   0.0621  0.0488   0.0741   0.0457  0.0590   0.0433   0.0412  0.1274
##          J2L2     J3L1    J3L2
##       -0.0029  -0.0028  0.0047
## s.e.   0.1243   0.0337  0.0335
## 
## sigma^2 estimated as 6914:  log likelihood = -3321.29,  aic = 6666.59
## 
## $degrees_of_freedom
## [1] 557
## 
## $ttable
##      Estimate     SE  t.value p.value
## ar1   -0.1167 0.0621  -1.8780  0.0609
## ar2    0.8114 0.0488  16.6144  0.0000
## ma1   -0.4463 0.0741  -6.0209  0.0000
## ma2   -0.9757 0.0457 -21.3353  0.0000
## ma3    0.4221 0.0590   7.1540  0.0000
## sar1  -0.6771 0.0433 -15.6260  0.0000
## sar2  -0.2916 0.0412  -7.0873  0.0000
## J2L1   0.3528 0.1274   2.7701  0.0058
## J2L2  -0.0029 0.1243  -0.0233  0.9814
## J3L1  -0.0028 0.0337  -0.0840  0.9331
## J3L2   0.0047 0.0335   0.1391  0.8894
## 
## $AIC
## [1] 11.73695
## 
## $AICc
## [1] 11.73779
## 
## $BIC
## [1] 11.82869
## 
## model 13
## 
## Series: data.1$J1 
## Regression with ARIMA(1,0,1)(0,1,1)[7] errors 
## 
## Coefficients:
##          ar1      ma1     sma1   drift    J2L1     J2L2     J3L1   J3L2
##       0.8234  -0.4002  -0.8220  1.9402  0.3404  -0.0569  -0.0028  0.003
## s.e.  0.0390   0.0670   0.0249  0.3212  0.1251   0.1191   0.0354  0.035
## 
## sigma^2 = 6300:  log likelihood = -3296.21
## AIC=6610.41   AICc=6610.74   BIC=6649.51
## model 14n
## Series: data.1$J1 
## Regression with ARIMA(1,0,1)(0,1,1)[7] errors 
## 
## Coefficients:
##          ar1      ma1     sma1   drift    xreg
##       0.8220  -0.4012  -0.8206  1.9064  0.3315
## s.e.  0.0394   0.0669   0.0248  0.3123  0.1228
## 
## sigma^2 = 6270:  log likelihood = -3296.32
## AIC=6604.65   AICc=6604.8   BIC=6630.71

\[ (1-0.822B)(1+0.82B^7)\nabla^{7}J_{1,t} = 1.906 +(1-0.401B)w_{1,t} + J_{2,t-1} \]

J1 model evaluation

hist(model1.14$residuals, main = "histogram of residuals", freq = F, breaks = 35)
lines(-400:400,dnorm(-400:400, mean = mean(model1.14$residuals), sd = sd(model1.14$residuals)), col = "red")

plot(model1.14$residuals, main = "residual plot", type = "l")

adf.test(model1.14$residuals, alternative = "stationary")
## Warning in adf.test(model1.14$residuals, alternative = "stationary"): p-value
## smaller than printed p-value
acf(model1.14$residuals, main = "model 14 residual ACF")
pacf(model1.14$residuals, main = "model 14 residual PACF")


model1.forecast1 = forecast(cars.day.fit1, h = 30)
plot(model1.forecast1, xlim = c(75,90))
lines(cars.day$Junction1)


model1.forecast2 = forecast(model1.14, h = 30, xreg = data.1.full$J2L1[577:606])
plot(model1.forecast2, xlim = c(560,610))
lines(data.1.full$J1)

rmse1.1 = rmse(model1.forecast1$mean, data.1.full$J1[577:606])
# 114.3474

rmse1.2 = rmse(model1.forecast2$mean, data.1.full$J1[577:606])
# 109.6458

# mean(data.1.full$J1[577:606])
# 1761.767
## 
##  Augmented Dickey-Fuller Test
## 
## data:  model1.14$residuals
## Dickey-Fuller = -8.2807, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary

Junction 2

#quiet(astsa::sarima(cars.day.train$Junction2, p = 3, d = 1, q = 3))

# ARIMA(2,0,1)(0,1,1)[7] with drift
# AIC 5435.15

data.2 = data.frame(J2 = cars.day.train$Junction2[-c(1,2)],
                    J1L1 = cars.day.train$Junction1[-c(1,578)],
                    J1L2 = cars.day.train$Junction1[-c(577,578)],
                    J3L1 = cars.day.train$Junction3[-c(1,578)],
                    J3L2 = cars.day.train$Junction3[-c(577,578)])

data.2.full = data.frame(J2 = cars.day$Junction2[-c(1,2)],
                         J1L1 = cars.day$Junction1[-c(1,608)],
                         J1L2 = cars.day$Junction1[-c(607,608)],
                         J3L1 = cars.day$Junction3[-c(1,608)],
                         J3L2 = cars.day$Junction3[-c(607,608)])



# (cars.day.fit2 = auto.arima(data.2$J2,
#                            max.p = 3,
#                         max.q = 3,
#                         max.P = 3,
#                         max.Q = 3,
#                         max.order = 5,
#                         max.d = 1,
#                         max.D = 1,
#                         start.p = 2,
#                         start.q = 1,
#                         start.P = 1,
#                         start.Q = 1,
#                         nmodels = 100,
#                         approximation = T, 
#                         seasonal = T))



(model2.1 = Arima(data.2$J2, 
          order=c(2, 0, 1), 
          seasonal=list(order=c(0, 1, 1), period=7), 
          include.drift=T,
          xreg = as.matrix(data.2[,-c(1)])))
#AIC 5415.31

(model2.2 = Arima(data.2$J2, 
          order=c(2, 0, 1), 
          seasonal=list(order=c(0, 1, 1), period=7), 
          include.drift=T,
          xreg = as.matrix(data.2[,-c(1,3)])))
#AIC 5413.33

(model2.3 = Arima(data.2$J2, 
          order=c(2, 0, 1), 
          seasonal=list(order=c(0, 1, 1), period=7), 
          include.drift=T,
          xreg = as.matrix(data.2[,-c(1,3,5)])))
#AIC 5411.81
## Series: data.2$J2 
## Regression with ARIMA(2,0,1)(0,1,1)[7] errors 
## 
## Coefficients:
##          ar1      ar2      ma1     sma1   drift    J1L1    J1L2    J3L1    J3L2
##       1.3769  -0.3915  -0.8586  -0.7987  0.6454  0.0301  0.0019  0.0260  0.0080
## s.e.  0.0598   0.0558   0.0365   0.0264  0.3212  0.0157  0.0150  0.0122  0.0122
## 
## sigma^2 = 771.7:  log likelihood = -2697.66
## AIC=5415.31   AICc=5415.71   BIC=5458.75
## Series: data.2$J2 
## Regression with ARIMA(2,0,1)(0,1,1)[7] errors 
## 
## Coefficients:
##          ar1      ar2      ma1     sma1   drift    J1L1    J3L1    J3L2
##       1.3771  -0.3917  -0.8580  -0.7990  0.6500  0.0299  0.0260  0.0083
## s.e.  0.0600   0.0560   0.0363   0.0263  0.3193  0.0157  0.0122  0.0119
## 
## sigma^2 = 770.3:  log likelihood = -2697.66
## AIC=5413.33   AICc=5413.65   BIC=5452.42
## Series: data.2$J2 
## Regression with ARIMA(2,0,1)(0,1,1)[7] errors 
## 
## Coefficients:
##          ar1      ar2      ma1     sma1   drift    J1L1    J3L1
##       1.3760  -0.3912  -0.8552  -0.7990  0.6554  0.0292  0.0278
## s.e.  0.0611   0.0569   0.0373   0.0263  0.3151  0.0157  0.0119
## 
## sigma^2 = 769.6:  log likelihood = -2697.9
## AIC=5411.81   AICc=5412.06   BIC=5446.56

Junction 2 forecasts

# ARIMA(2,0,1)(0,1,1)[7] with drift 
model2.forecast1 = forecast(cars.day.fit2, h = 30)
plot(model2.forecast1, xlim = c(75,90))
lines(cars.day$Junction2)
(rmse2.1 = rmse(model2.forecast1$mean, data.2.full[577:606,1]))
# 43.50
# AIC 5435.15

model2.forecast2 = forecast(model2.1, h = 30, 
                            xreg = as.matrix(data.2.full[577:606,-c(1)]))
plot(model2.forecast2, xlim = c(560,610))
lines(data.2.full$J2)
(rmse2.2 = rmse(model2.forecast2$mean, data.2.full[577:606,1]))
# 43.36

model2.forecast3 = forecast(model2.3, h = 30, 
                            xreg = as.matrix(data.2.full[577:606,-c(1,3,5)]))
plot(model2.forecast3, xlim = c(560,610))
lines(data.2.full$J2)
(rmse2.2 = rmse(model2.forecast3$mean, data.2.full[577:606,1]))
#42.88
## [1] 43.50186
## [1] 43.36398
## [1] 42.88388

Junction 3

#quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d = 1, q = 1))

#ARIMA(1,1,1)(0,0,1)[7] 
# AIC 6878.07

data.3 = data.frame(J3 = cars.day.train$Junction3[-c(1,2)],
                    J1L1 = cars.day.train$Junction1[-c(1,578)],
                    J1L2 = cars.day.train$Junction1[-c(577,578)],
                    J2L1 = cars.day.train$Junction2[-c(1,578)],
                    J2L2 = cars.day.train$Junction2[-c(577,578)])

data.3.full = data.frame(J3 = cars.day$Junction3[-c(1,2)],
                         J1L1 = cars.day$Junction1[-c(1,608)],
                         J1L2 = cars.day$Junction1[-c(607,608)],
                         J2L1 = cars.day$Junction2[-c(1,608)],
                         J2L2 = cars.day$Junction2[-c(607,608)])


# cars.day.fit3
# ARIMA(1,1,1)(0,0,1)[7] 
# AIC=6878.07 


(model3.1 = Arima(data.3$J3, 
          order=c(1, 1, 1), 
          seasonal=list(order=c(0, 0, 1), period=7), 
          include.drift=T,
          xreg = as.matrix(data.3[,-c(1)]),
          method = "ML"))
#AIC 6857.59

(model3.2 = Arima(data.3$J3, 
          order=c(1, 1, 1), 
          seasonal=list(order=c(0, 0, 1), period=7), 
          include.drift=T,
          xreg = as.matrix(data.3[,-c(1,3,5)]),
          method = "ML"))
#AIC 6854.43

(model3.3 = Arima(data.3$J3, 
          order=c(1, 1, 1), 
          seasonal=list(order=c(0, 0, 1), period=7), 
          include.drift=T,
          xreg = as.matrix(data.3[,-c(1,2,3,5)]),
          method = "ML"))
#AIC 6853.18

(model3.4 = Arima(data.3$J3, 
          order=c(1, 1, 1), 
          seasonal=list(order=c(1, 0, 1), period=7), 
          include.drift=T,
          #xreg = as.matrix(data.3[,-c(1,2,3,5)]),
          method = "ML"))
#AIC 6842.05
## Series: data.3$J3 
## Regression with ARIMA(1,1,1)(0,0,1)[7] errors 
## 
## Coefficients:
##          ar1      ma1    sma1   drift     J1L1    J1L2    J2L1     J2L2
##       0.7154  -1.0000  0.0656  0.4232  -0.0304  0.0233  0.1806  -0.1047
## s.e.  0.0298   0.0063  0.0419  0.1056   0.0386  0.0387  0.1183   0.1179
## 
## sigma^2 = 8632:  log likelihood = -3419.79
## AIC=6857.59   AICc=6857.91   BIC=6896.78
## Series: data.3$J3 
## Regression with ARIMA(1,1,1)(0,0,1)[7] errors 
## 
## Coefficients:
##          ar1      ma1    sma1   drift     J1L1    J2L1
##       0.7156  -1.0000  0.0671  0.4167  -0.0330  0.1768
## s.e.  0.0298   0.0069  0.0417  0.0969   0.0381  0.1157
## 
## sigma^2 = 8615:  log likelihood = -3420.22
## AIC=6854.43   AICc=6854.63   BIC=6884.91
## Series: data.3$J3 
## Regression with ARIMA(1,1,1)(0,0,1)[7] errors 
## 
## Coefficients:
##          ar1      ma1    sma1   drift    xreg
##       0.7133  -1.0000  0.0683  0.3975  0.0932
## s.e.  0.0297   0.0064  0.0417  0.0939  0.0639
## 
## sigma^2 = 8611:  log likelihood = -3420.59
## AIC=6853.18   AICc=6853.33   BIC=6879.31
## Series: data.3$J3 
## ARIMA(1,1,1)(1,0,1)[7] with drift 
## 
## Coefficients:
##          ar1      ma1    sar1     sma1   drift
##       0.7321  -0.9999  0.9917  -0.9711  0.4593
## s.e.  0.0286   0.0037  0.0092   0.0180  0.0991
## 
## sigma^2 = 8429:  log likelihood = -3415.03
## AIC=6842.05   AICc=6842.2   BIC=6868.18

J3 forecasting

model3.forecast1 = forecast(cars.day.fit3, h = 30)
plot(model3.forecast1, xlim = c(75,90))
lines(cars.day$Junction3)
(rmse3.1 = rmse(model3.forecast1$mean, data.3.full[577:606,1]))
# 124.8629


model3.forecast2 = forecast(model3.1, h = 30, 
                            xreg = as.matrix(data.3.full[577:606,-c(1)]))
plot(model3.forecast2, xlim = c(560,610))
lines(data.3.full$J3)
(rmse3.2 = rmse(model3.forecast2$mean, data.3.full[577:606,1]))
# 110.0729

model3.forecast3 = forecast(model3.4, h = 30)
plot(model3.forecast3, xlim = c(560,610))
lines(data.3.full$J3)
(rmse3.2 = rmse(model3.forecast3$mean, data.3.full[577:606,1]))
# 114.9068
## [1] 124.8629
## [1] 110.0729
## [1] 114.9068

look at residuals for better fit

#Hi

TBATS

fit <- tbats(data.1$J1, seasonal.periods = c(7))
fore = forecast(fit, h = 30)
plot(fore, xlim = c(75,90));
lines(ts(data.1.full$J1,frequency = 7))
# cat("model AIC is", fit$AIC) 
# AIC 8661.047
rmse.tbats = rmse(data.1.full$J1[577:606], fore$mean)
#148.5043

Log fits

J1 log

(model1.log.1 = auto.arima(log(data.1$J1),
                           max.p = 5,
                        max.q = 5,
                        max.P = 5,
                        max.Q = 5,
                        max.order = 10,
                        max.d = 3,
                        max.D = 3,
                        start.p = 1,
                        start.q = 1,
                        start.P = 1,
                        start.Q = 1,
                        nmodels = 100,
                        approximation = T))

model1.log.forecast1 = forecast(model1.log.1, h = 30)
plot(model3.forecast3)
#, xlim = c(560,610)

lines(log(data.1.full$J1))
(rmse1.log.1 = rmse(model1.log.forecast1$mean, data.1.full[577:606,1]))
#1790.762
## Series: log(data.1$J1) 
## ARIMA(5,1,1) with drift 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ar5      ma1   drift
##       -0.0768  -0.5805  -0.3777  -0.3406  -0.6302  -0.4245  0.0024
## s.e.   0.0371   0.0292   0.0362   0.0295   0.0341   0.0349  0.0010
## 
## sigma^2 = 0.01619:  log likelihood = 370.64
## AIC=-725.27   AICc=-725.02   BIC=-690.44
## [1] 1790.762

cat("1\n\n\n\n")

auto.arima(cars.day.train$Junction1,
                           max.p = 5,
                           max.q = 5,
                           max.P = 5,
                           max.Q = 5,
                           max.order = 10,
                           max.d = 3,
                           max.D = 3,
                           start.p = 1,
                           start.q = 1,
                           start.P = 1,
                           start.Q = 1,
                           nmodels = 100,
                           approximation = T)

cat("\n\n2\n\n\n\n")

auto.arima(cars.day.train$Junction1[-c(1,2)],
                           max.p = 5,
                           max.q = 5,
                           max.P = 5,
                           max.Q = 5,
                           max.order = 10,
                           max.d = 3,
                           max.D = 3,
                           start.p = 1,
                           start.q = 1,
                           start.P = 1,
                           start.Q = 1,
                           nmodels = 100,
                           approximation = T)
## 1
## 
## 
## 
## Series: cars.day.train$Junction1 
## ARIMA(1,0,1)(0,1,1)[7] with drift 
## 
## Coefficients:
##          ar1      ma1     sma1   drift
##       0.8128  -0.3473  -0.8211  2.1295
## s.e.  0.0386   0.0629   0.0243  0.3114
## 
## sigma^2 = 6325:  log likelihood = -3310.96
## AIC=6631.91   AICc=6632.02   BIC=6653.65
## 
## 
## 2
## 
## 
## 
## Series: cars.day.train$Junction1[-c(1, 2)] 
## ARIMA(3,1,3) 
## 
## Coefficients:
##          ar1     ar2      ar3      ma1      ma2     ma3
##       0.2874  0.0921  -0.5859  -0.6139  -0.6891  0.6527
## s.e.  0.0636  0.0648   0.0657   0.0503   0.0661  0.0519
## 
## sigma^2 = 24274:  log likelihood = -3717.2
## AIC=7448.39   AICc=7448.59   BIC=7478.87
my.mat1 = matrix(0, nrow = 2, ncol = 3)
colnames(my.mat1) = c("Model 1", "Model 2", "Model 3")
rownames(my.mat1) = c("AIC","RMSE")

my.mat1[1,1] = 6631.91
my.mat1[1,2] = 6604.65
my.mat1[1,3] = 6608.14
my.mat1[2,1] = 114.34
my.mat1[2,2] = 109.64
my.mat1[2,3] = 115.10

knitr::kable(my.mat1)


my.mat2 = matrix(0, nrow = 2, ncol = 3)
colnames(my.mat2) = c("Model 1", "Model 2", "Model 3")
rownames(my.mat2) = c("AIC","RMSE")

my.mat2[1,1] = 5435.15
my.mat2[1,2] = 5415.31
my.mat2[1,3] = 5411.81
my.mat2[2,1] = 43.50
my.mat2[2,2] = 43.36
my.mat2[2,3] = 42.88

knitr::kable(my.mat2)



my.mat3 = matrix(0, nrow = 2, ncol = 3)
colnames(my.mat3) = c("Model 1", "Model 2", "Model 3")
rownames(my.mat3) = c("AIC","RMSE")

my.mat3[1,1] = 6878.07
my.mat3[1,2] = 6857.59
my.mat3[1,3] = 6842.05
my.mat3[2,1] = 124.86
my.mat3[2,2] = 110.07
my.mat3[2,3] = 114.90

knitr::kable(my.mat3)
Model 1 Model 2 Model 3
AIC 6631.91 6604.65 6608.14
RMSE 114.34 109.64 115.10
Model 1 Model 2 Model 3
AIC 5435.15 5415.31 5411.81
RMSE 43.50 43.36 42.88
Model 1 Model 2 Model 3
AIC 6878.07 6857.59 6842.05
RMSE 124.86 110.07 114.90

note plot residuals for all models

LateX

Model 1

1

\[ (1-0.81B)(1+0.82B^7)(1-B^7)J_{1,t} = 2.12 + (1-0.34B)w_{1,t} \]

par(mfrow = c(2,2))

plot(cars.day.fit1$residuals, ylab = "residual", main = "residuals")
hist(cars.day.fit1$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(cars.day.fit1$residuals), sd = sd(cars.day.fit1$residuals)), col = "red")
acf(cars.day.fit1$residuals, lag.max = 15, main = "residual ACF")
pacf(cars.day.fit1$residuals, lag.max = 15, main = "residual PACF")

2

\[ (1-0.822B)(1+0.82B^7)(1-B^7)J_{1,t} = 1.906 +(1-0.401B)w_{1,t} + 0.33J_{2,t-1} \]

par(mfrow = c(2,2))

plot(model1.14$residuals, ylab = "residual", main = "residuals")
hist(model1.14$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model1.14$residuals), sd = sd(model1.14$residuals)), col = "red")
acf(model1.14$residuals, lag.max = 15, main = "residual ACF")
pacf(model1.14$residuals, lag.max = 15, main = "residual PACF")

\[ 1=2-1 \]

Model 2

\[ (1-1.41B+0.43B^2)(1-B^7)J_{2,t} = 0.72 + (1-0.84B)(1-0.79B^7)w_{2,t} \]

par(mfrow = c(2,2))

plot(cars.day.fit2$residuals, ylab = "residual", main = "residuals")
hist(cars.day.fit2$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(cars.day.fit2$residuals), sd = sd(cars.day.fit2$residuals)), col = "red")
acf(cars.day.fit2$residuals, lag.max = 15, main = "residual ACF")
pacf(cars.day.fit2$residuals, lag.max = 15, main = "residual PACF")

\[ (1-1.37B+0.39B^2)(1-B^7)J_{2,t} = 0.64 + (1-0.85B)(1-0.79B^7)w_{2,t} + 0.030J_{1,t-1} + 0.0019J_{1,t-2} + 0.026J_{3,t-1} + 0.0080J_{3,t-2} \]

par(mfrow = c(2,2))

plot(model2.1$residuals, ylab = "residual", main = "residuals")
hist(model2.1$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model2.1$residuals), sd = sd(model2.1$residuals)), col = "red")
acf(model2.1$residuals, lag.max = 15, main = "residual ACF")
pacf(model2.1$residuals, lag.max = 15, main = "residual PACF")

\[ (1-1.37B+0.39B^2)(1-B^7)J_{2,t} = 0.65 + (1-0.85B)(1-0.79B^7)w_{3,t} +0.029J_{1,t-1} 0.027J_{3,t-1} \]

par(mfrow = c(2,2))

plot(model2.3$residuals, ylab = "residual", main = "residuals")
hist(model2.3$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model2.3$residuals), sd = sd(model2.3$residuals)), col = "red")
acf(model2.3$residuals, lag.max = 15, main = "residual ACF")
pacf(model2.3$residuals, lag.max = 15, main = "residual PACF")

Model 3

\[ (1-0.69B)\nabla J_{3,t} = (1-0.97B)(1+0.067B^7)w_{3,t} \]

par(mfrow = c(2,2))

plot(cars.day.fit3$residuals, ylab = "residual", main = "residuals")
hist(cars.day.fit3$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(cars.day.fit3$residuals), sd = sd(cars.day.fit3$residuals)), col = "red")
acf(cars.day.fit3$residuals, lag.max = 15, main = "residual ACF")
pacf(cars.day.fit3$residuals, lag.max = 15, main = "residual PACF")

\[ (1-0.71B)\nabla J_{3,t} = 0.42+(1-1B)(1+0.065B^7)w_{3,t} -0.030J_{1,t-1} + 0.023J_{1,t-2} + 0.18J_{2,t-1} -0.10J_{2,t-2} \]

par(mfrow = c(2,2))

plot(model3.1$residuals, ylab = "residual", main = "residuals")
hist(model3.1$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model3.1$residuals), sd = sd(model3.1$residuals)), col = "red")
acf(model3.1$residuals, lag.max = 15, main = "residual ACF")
pacf(model3.1$residuals, lag.max = 15, main = "residual PACF")

\[ (1-0.73B)(1-0.99B^7)\nabla J_{3,t} = 0.45 + (1-B)(1-0.97B^7)w_{3,t} \]

par(mfrow = c(2,2))

plot(model3.4$residuals, ylab = "residual", main = "residuals")
hist(model3.4$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model3.4$residuals), sd = sd(model3.4$residuals)), col = "red")
acf(model3.4$residuals, lag.max = 15, main = "residual ACF")
pacf(model3.4$residuals, lag.max = 15, main = "residual PACF")

Relationships between vars in the data.numbers

J1

par(mfrow = c(2,2))
plot(data.1.full$J2L1,data.1.full$J1, main = "J1 | J2(t-1)", xlab = "J2(t-1)", ylab = "J1")
abline(v = 450, col = "red", lty = 2)
plot(data.1.full$J2L2,data.1.full$J1, main = "J1 | J2(t-2)", xlab = "J2(t-2)", ylab = "J1")
abline(v = 450, col = "red", lty = 2)
plot(data.1.full$J3L1,data.1.full$J1, main = "J1 | J3(t-1)", xlab = "J3(t-1)", ylab = "J1")
plot(data.1.full$J3L2,data.1.full$J1, main = "J1 | J3(t-2)", xlab = "J3(t-2)", ylab = "J1")

data.1.full2 = data.1.full

data.1.full2$J2L1L450 = ifelse(data.1.full$J2L1<450,data.1.full$J2L1,0)
data.1.full2$J2L1G450 = ifelse(data.1.full$J2L1>=450,1,0)
#data.1.full2$J2L1L4502 = ifelse(data.1.full$J2L1<450,data.1.full$J2L1,0)


data.1.2 = data.1.full2[-c(577:606),]

(model1.15 = Arima(data.1$J1, 
          order=c(1, 0, 1), 
          seasonal=list(order=c(0, 1, 1), period=7), 
          include.drift=T,
          xreg = as.matrix(data.1.2[,-c(1,2,3,4,5)])))

#AIC 6608.14

par(mfrow = c(2,2))
plot(model1.15$residuals, ylab = "residual", main = "residuals")
hist(model1.15$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model1.15$residuals), sd = sd(model1.15$residuals)), col = "red")
acf(model1.15$residuals, lag.max = 15, main = "residual ACF")
pacf(model1.15$residuals, lag.max = 15, main = "residual PACF")


par(mfrow = c(1,1))
model1.forecast3 = forecast(model1.15, h = 30, xreg = as.matrix(data.1.full2[577:606,-c(1:5)]))
plot(model1.forecast2, xlim = c(560,610))
lines(data.1.full$J1)

rmse1.3 = rmse(model1.forecast3$mean, data.1.full$J1[577:606])
#115.1091
## Series: data.1$J1 
## Regression with ARIMA(1,0,1)(0,1,1)[7] errors 
## 
## Coefficients:
##          ar1      ma1     sma1   drift  J2L1L450  J2L1G450
##       0.8067  -0.3697  -0.8134  1.9971    0.3223  143.3700
## s.e.  0.0413   0.0662   0.0260  0.3110    0.1350   67.5697
## 
## sigma^2 = 6300:  log likelihood = -3297.07
## AIC=6608.14   AICc=6608.34   BIC=6638.55

\[ (1-0.80B)(1-B^7)J_{1,t} = 1.99+(1-0.36B)(1-0.81B^7)+ 0.32J_{2,t-1}(J_{2,t-1}<450) + 143.37(J_{2,t-1}\ge450) \]

test \(\mathbb{I}\)
\[ \mathbb{I} \]

J2

par(mfrow = c(2,2))
plot(data.2.full$J1L1,data.2.full$J2, main = "J2 | J1(t-1)", xlab = "J1(t-1)", ylab = "J2")
plot(data.2.full$J1L2,data.2.full$J2, main = "J2 | J1(t-2)", xlab = "J1(t-2)", ylab = "J2")
plot(data.2.full$J3L1,data.2.full$J2, main = "J2 | J3(t-1)", xlab = "J3(t-1)", ylab = "J2")
plot(data.2.full$J3L2,data.2.full$J2, main = "J2 | J3(t-2)", xlab = "J3(t-2)", ylab = "J2")

J3

par(mfrow = c(2,2))
plot(data.3.full$J1L1,data.3.full$J3, main = "J3 | J1(t-1)", xlab = "J1(t-1)", ylab = "J3")
plot(data.3.full$J1L2,data.3.full$J3, main = "J3 | J1(t-2)", xlab = "J1(t-2)", ylab = "J3")
plot(data.3.full$J2L1,data.3.full$J3, main = "J3 | J2(t-1)", xlab = "J2(t-1)", ylab = "J3")
plot(data.3.full$J2L2,data.3.full$J3, main = "J3 | J2(t-2)", xlab = "J2(t-2)", ylab = "J3")

residulal correlations

# J1 model1.14
# J2 model2.3
# J3 model3.4

ggCcf(model1.14$residuals, model2.3$residuals)+ggtitle("CCF Plot for Junctions 1 & 2")
ggCcf(model1.14$residuals, model3.4$residuals)+ggtitle("CCF Plot for Junctions 1 & 3")
ggCcf(model2.3$residuals, model3.4$residuals)+ggtitle("CCF Plot for Junctions 2 & 3")

\[ (1-0.81B)(1-B^7)J_{1,t} = (1-0.34B)(1-0.82B^7)w_{1,t} \]

\[ (1-0.82B)(1-B^7)J_{1,t} = 1.90+ (1-0.40B)(1-0.82B^7)w_{1,t} +0.33J_{2,t-1} \]